Multinomial Logistic Regression
Advanced Categorical Data Analysis
This document serves both learning and practical purposes. It is designed for educational use, aiming to enhance statistical analysis skills and provide clear, organized notes for future reference.
1 Introduction
Multinomial Logistic Regression is an extension of binary logistic regression that allows for the modeling of outcome variables with more than two categories. This type of regression is used when the outcome categories are nominal, meaning they do not have a natural order. The dependent variable has multiple categories (more than two), and one of the categories is chosen as the reference category. The model compares each of the other categories to this reference category.
1.1 Model Structure
- The model estimates the log odds of each category relative to the reference category.
- For a categorical outcome with \(k\) categories, the model fits \(k-1\) equations.
- Each equation models the log odds of a specific category relative to the reference category.
1.2 Model Formulation
- Log Odds
- The log odds of being in category \(j\) relative to the reference category \(0\) is modeled as:
\[ \log \left( \frac{P(Y=j|X)}{P(Y=0|X)} \right) = \alpha_j + \beta_{j1}X_1 + \beta_{j2}X_2 + \ldots + \beta_{jp}X_p \]
\(\alpha_j\) is the intercept for category \(j\), and \(\beta_{jk}\) are the coefficients for the predictor variables \(X_k\).
- Probability
- The probability of being in category \(j\) is derived from the log odds:
\[ P(Y=j|X) = \frac{e^{\alpha_j + \beta_{j1}X_1 + \beta_{j2}X_2 + \ldots + \beta_{jp}X_p}}{1 + \sum_{l=1}^{k-1} e^{\alpha_l + \beta_{l1}X_1 + \beta_{l2}X_2 + \ldots + \beta_{lp}X_p}} \]
- The probability of being in the reference category \(0\) is:
\[ P(Y=0|X) = \frac{1}{1 + \sum_{l=1}^{k-1} e^{\alpha_l + \beta_{l1}X_1 + \beta_{l2}X_2 + \ldots + \beta_{lp}X_p}} \]
1.3 Estimation
- The parameters (\(\alpha_j\) and \(\beta_{jk}\)) are estimated using maximum likelihood estimation (MLE).
- The likelihood function is constructed based on the observed data and the probabilities of each category.
- Iterative algorithms, such as Newton-Raphson or Fisher scoring, are used to find the parameter estimates that maximize the likelihood function.
1.4 Interpretation
-
Coefficients (\(\beta\))
- The coefficients represent the change in the log odds of being in a category relative to the reference category for a one-unit change in the predictor variable.
- Positive coefficients indicate higher odds of being in the category compared to the reference, while negative coefficients indicate lower odds.
Odds Ratios (OR) is obtained by exponentiating the coefficients:
\[ OR = e^{\beta} \]
An OR greater than 1 indicates increased odds of being in the category compared to the reference, and an OR less than 1 indicates decreased odds.
1.5 Model Fit and Diagnostics
- Likelihood Ratio Test: Used to compare nested models and assess the significance of adding or removing predictors.
- Wald Test: Used to test the significance of individual coefficients.
2 Setting Up the Environment
Load the necessary libraries to ensure the R environment is equipped with the tools and functions required for efficient and effective analysis.
3 Practical 1: Mammography Experience Dataset
In this exercise, the objectives are to understand and apply multinomial logistic regression using a dataset related to mammography experiences. The focus is on analyzing the mammography experiences dataset to understand factors influencing whether women have had a mammogram recently, a long time ago, or never.
3.1 Dataset Overview
The dataset used is focused on mammography experiences which consists of responses from a survey on mammography experiences and perceptions. The dataset includes the following variables:
-
me (Mammography Experience): This is the outcome variable and is categorical with three levels:
-
0: Never had a mammogram -
1: Had a mammogram within a year -
2: Had a mammogram over a year ago
-
-
sympt (Symptoms Perception): This variable captures the perception that a mammogram is not needed unless symptoms are present.
-
1: Strongly Agree -
2: Agree -
3: Disagree -
4: Strongly Disagree
-
pb (Perceived Benefit of Mammography): This variable measures the perceived benefit of mammography on a scale. The exact scale is not specified but typically ranges from lower to higher perceived benefits.
-
hist (Family History of Breast Cancer): This binary variable indicates whether the respondent has a family history (mother or sister) of breast cancer.
-
0: No -
1: Yes
-
-
bse (Breast Self-Examination Training): This binary variable indicates whether the respondent has been taught how to perform breast self-examinations.
-
0: No -
1: Yes
-
-
detc (Detection Likelihood): This variable measures the respondent’s perception of the likelihood that a mammogram could detect a new case of breast cancer. It is coded as:
-
1: Not likely -
2: Somewhat likely -
3: Very likely
-
3.2 Data Reading and Exploration
Rows: 412
Columns: 7
$ obs <dbl> 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 1…
$ me <dbl> 0, 0, 0, 1, 2, 0, 2, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 2, 0…
$ sympt <dbl> 3, 2, 3, 3, 4, 3, 4, 4, 2, 4, 4, 3, 4, 1, 2, 4, 3, 4, 3, 3, 3, 4…
$ pb <dbl> 7, 11, 8, 11, 7, 7, 6, 6, 6, 6, 8, 6, 6, 5, 8, 11, 6, 5, 10, 10,…
$ hist <dbl> 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0…
$ bse <dbl> 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 0, 1, 1, 0, 1, 1, 1, 1, 1, 1, 1…
$ detc <dbl> 2, 3, 3, 3, 3, 3, 2, 3, 3, 3, 2, 3, 3, 3, 2, 3, 3, 3, 3, 3, 3, 2…summary(dat.m1) obs me sympt pb
Min. : 1.0 Min. :0.0000 Min. :1.000 Min. : 5.000
1st Qu.:103.8 1st Qu.:0.0000 1st Qu.:2.000 1st Qu.: 6.000
Median :206.5 Median :0.0000 Median :3.000 Median : 7.000
Mean :206.5 Mean :0.6117 Mean :2.966 Mean : 7.561
3rd Qu.:309.2 3rd Qu.:1.0000 3rd Qu.:4.000 3rd Qu.: 9.000
Max. :412.0 Max. :2.0000 Max. :4.000 Max. :17.000
hist bse detc
Min. :0.0000 Min. :0.0000 Min. :1.000
1st Qu.:0.0000 1st Qu.:1.0000 1st Qu.:2.000
Median :0.0000 Median :1.0000 Median :3.000
Mean :0.1068 Mean :0.8689 Mean :2.658
3rd Qu.:0.0000 3rd Qu.:1.0000 3rd Qu.:3.000
Max. :1.0000 Max. :1.0000 Max. :3.000 3.3 Data Wrangling
The me and hist variables are treated as numerical despite being categorical variables. To ensure that the data is handled correctly during statistical analyses, these variables must be converted to factors so that each level is treated as a distinct category, and the analysis methods can correctly interpret these categories.
3.3.1 me Variable
# A tibble: 3 × 2
me n
<dbl> <int>
1 0 234
2 1 104
3 2 74# Copy dataset and keep the original dataset
dat.m2 <- dat.m1
# Converting `me` to a factor
dat.m2 <- dat.m2 %>%
mutate(me2 = factor(me, labels = c("never", "within.a.year", "over.a.year.ago")))
dat.m2 %>% count(me2)# A tibble: 3 × 2
me2 n
<fct> <int>
1 never 234
2 within.a.year 104
3 over.a.year.ago 74This code creates a new variable me2 by converting the numeric me variable to a factor. The labels argument assigns meaningful labels to the factor levels.
3.3.2 hist Variable
# A tibble: 2 × 2
hist n
<dbl> <int>
1 0 368
2 1 44# Converting `hist` to a Factor
dat.m2 <- dat.m2 %>% mutate(hist2 = factor(hist, labels = c("no", "yes")))
dat.m2 %>% count(hist2)# A tibble: 2 × 2
hist2 n
<fct> <int>
1 no 368
2 yes 44This code creates a new variable hist2 by converting the numeric hist variable to a factor, with no and yes as labels.
3.3.3 Table Summary
# Generate summary statistics for the dataset, grouped by the `me2` variable
dat.m2 %>%
select(-me, -hist) %>%
tbl_summary(by = me2, statistic = list(all_continuous() ~ "{mean} ({sd})"))| Characteristic | never, N = 2341 | within.a.year, N = 1041 | over.a.year.ago, N = 741 |
|---|---|---|---|
| obs | 203 (121) | 212 (113) | 211 (122) |
| sympt | |||
| 1 | 33 (14%) | 2 (1.9%) | 5 (6.8%) |
| 2 | 62 (26%) | 4 (3.8%) | 7 (9.5%) |
| 3 | 85 (36%) | 43 (41%) | 32 (43%) |
| 4 | 54 (23%) | 55 (53%) | 30 (41%) |
| pb | 8.06 (2.19) | 6.70 (1.66) | 7.19 (1.91) |
| bse | 190 (81%) | 99 (95%) | 69 (93%) |
| detc | |||
| 1 | 13 (5.6%) | 1 (1.0%) | 4 (5.4%) |
| 2 | 77 (33%) | 12 (12%) | 16 (22%) |
| 3 | 144 (62%) | 91 (88%) | 54 (73%) |
| hist2 | 14 (6.0%) | 19 (18%) | 11 (15%) |
| 1 Mean (SD); n (%) | |||
3.4 Estimation
The VGAM package is used to run multinomial logistic regression, which is suitable for a dependent variable with more than two categories. The objctive is to use the VGAM package to fit a multinomial logistic regression model and ensure the reference category for the dependent variable me2 is “never”.
3.4.1 Confirm the Order of Factor Levels
Check the current order of the factor levels for me2. It is crucial because the reference level in regression models affects the interpretation of coefficients.
levels(dat.m2$me2)[1] "never" "within.a.year" "over.a.year.ago"The levels are currently in the order where “never” is the first level. For the model, we need “never” to be the last level to serve as the reference category.
3.4.2 Reorder the Factor Levels
# Create new dataset `dat.m3`
dat.m3 <- dat.m2
# Generate `me3`
dat.m3 <- dat.m3 %>%
mutate(me3 = fct_relevel(me2, c("over.a.year.ago", "within.a.year", "never")))Creates a new variable me3 by reordering the levels of me2 using fct_relevel so that “never” becomes the reference level.
table(dat.m3$me)
0 1 2
234 104 74 summary(dat.m3$me2) never within.a.year over.a.year.ago
234 104 74 summary(dat.m3$me3)over.a.year.ago within.a.year never
74 104 234 These functions confirm that the observations and levels are correctly adjusted after releveling.
3.4.3 Check the Order of Levels
levels(dat.m3$me2)[1] "never" "within.a.year" "over.a.year.ago"levels(dat.m3$me3)[1] "over.a.year.ago" "within.a.year" "never" This checks the new order of levels to ensure that me3 has “never” as the last level (reference category).
3.4.4 Multinomial Logistic Regression Model
3.4.4.1 Model with hist2 as Predictor
fitmlog1 <- vglm(me3 ~ hist2,
family = multinomial(),
data = dat.m3
)
summary(fitmlog1)
Call:
vglm(formula = me3 ~ hist2, family = multinomial(), data = dat.m3)
Coefficients:
Estimate Std. Error z value Pr(>|z|)
(Intercept):1 -1.2505 0.1429 -8.751 < 2e-16 ***
(Intercept):2 -0.9510 0.1277 -7.446 9.6e-14 ***
hist2yes:1 1.0093 0.4275 2.361 0.018225 *
hist2yes:2 1.2564 0.3747 3.353 0.000798 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Names of linear predictors: log(mu[,1]/mu[,3]), log(mu[,2]/mu[,3])
Residual deviance: 792.3399 on 820 degrees of freedom
Log-likelihood: -396.17 on 820 degrees of freedom
Number of Fisher scoring iterations: 4
No Hauck-Donner effect found in any of the estimates
Reference group is level 3 of the response-
vglm(me3 ~ hist2, family = multinomial(), data = dat.m1): Fits a model whereme3(mammography experience) is the dependent variable, andhist2(family history of breast cancer) is the predictor. -
family = multinomial(): Specifies that the model is a multinomial logistic regression. -
data = dat.m3: Specifies the dataset to be used.
Summary of Output
-
Intercepts
-
(Intercept):1: The log-odds of being in the “over a year ago” category relative to “never” whenhist2is the reference level (no family history). -
(Intercept):2: The log-odds of being in the “within a year” category relative to “never” whenhist2is the reference level (no family history).
-
-
Predictor (
hist2yes)-
hist2yes:1: The log-odds change of being in the “over a year ago” category relative to “never” for participants with a family history of breast cancer compared to those without. -
hist2yes:2: The log-odds change of being in the “within a year” category relative to “never” for participants with a family history of breast cancer compared to those without.
-
-
Statistical Significance
-
(Intercept):1and(Intercept):2are highly significant (p < 0.001), indicating that the intercepts are significantly different from zero. -
hist2yes:1is significant at the 0.05 level (*), indicating that the log-odds change for the “over a year ago” category relative to “never” is significant for those with a family history of breast cancer. -
hist2yes:2is highly significant (p < 0.001), indicating that the log-odds change for the “within a year” category relative to “never” is significant for those with a family history of breast cancer.
-
-
Interpretation
- The model indicates that having a family history of breast cancer (hist2 = “yes”) significantly affects the likelihood of having had a mammogram “within a year” or “over a year ago” compared to “never”.
- Significant positive coefficients for
hist2yessuggest that participants with a family history of breast cancer are more likely to have had a mammogram within a year or over a year ago compared to those without a family history. - The model fit statistics (residual deviance, log-likelihood) suggest that the model adequately fits the data.
3.4.4.2 Model with detc as Predictor
# Convert `detc` to a factor
dat.m3 <- dat.m3 %>% mutate(detc2 = factor(detc))
# Check the distribution of `detc2` across the levels of `me3`
table(dat.m3$me3, dat.m3$detc2)
1 2 3
over.a.year.ago 4 16 54
within.a.year 1 12 91
never 13 77 144The table shows the count of observations for each combination of me2 and detc2.
# Fitting the model with `detc2` as predictor
fitmlog2 <- vglm(me3 ~ detc2,
family = multinomial(),
data = dat.m3
)
summary(fitmlog2)
Call:
vglm(formula = me3 ~ detc2, family = multinomial(), data = dat.m3)
Coefficients:
Estimate Std. Error z value Pr(>|z|)
(Intercept):1 -1.1787 0.5718 -2.061 0.0393 *
(Intercept):2 -2.5649 1.0377 NA NA
detc22:1 -0.3926 0.6344 -0.619 0.5360
detc22:2 0.7061 1.0832 0.652 0.5145
detc23:1 0.1978 0.5936 0.333 0.7389
detc23:2 2.1060 1.0463 2.013 0.0441 *
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Names of linear predictors: log(mu[,1]/mu[,3]), log(mu[,2]/mu[,3])
Residual deviance: 778.4011 on 818 degrees of freedom
Log-likelihood: -389.2005 on 818 degrees of freedom
Number of Fisher scoring iterations: 5
Warning: Hauck-Donner effect detected in the following estimate(s):
'(Intercept):2'
Reference group is level 3 of the responseSummary of Output
-
Intercepts
-
(Intercept):1: The log-odds of being in the “over a year ago” category relative to “never” whendetc2is at its reference level (usually the first level). -
(Intercept):2: The log-odds of being in the “within a year” category relative to “never” whendetc2is at its reference level.
-
-
Predictor (
detc2)-
detc22:1: The log-odds change of being in the “over a year ago” category relative to “never” for participants with the second level ofdetc2compared to the reference level. -
detc22:2: The log-odds change of being in the “within a year” category relative to “never” for participants with the second level ofdetc2compared to the reference level. -
detc23:1: The log-odds change of being in the “over a year ago” category relative to “never” for participants with the third level ofdetc2compared to the reference level. -
detc23:2: The log-odds change of being in the “within a year” category relative to “never” for participants with the third level ofdetc2compared to the reference level.
-
-
**Statistical Significance*
-
(Intercept):1is significant at the 0.05 level (*), indicating that the intercept for the “over a year ago” category relative to “never” is significantly different from zero. -
detc23:2is significant at the 0.05 level (*), indicating that the log-odds change for the “within a year” category relative to “never” is significant for those at the third level ofdetc2.
-
-
Interpretation
- The model indicates that the perceived likelihood of detection (
detc2) has an effect on the likelihood of having had a mammogram “within a year” compared to “never”. - Significant positive coefficient for
detc23:2suggests that participants with the highest level of perceived likelihood of detection are more likely to have had a mammogram within a year compared to those with the lowest level of perceived likelihood. - The model fit statistics (residual deviance, log-likelihood) suggest that the model adequately fits the data.
- The model indicates that the perceived likelihood of detection (
3.5 Inferences
Inferences refer to the process of drawing conclusions about a population based on the analysis of sample data. In multinomial logistic regression model, inferences involve estimating the relationships between the predictor variables and the outcome variable and determining the statistical significance and confidence of these estimates. The process of inferences involves:
- Calculate the 95% Confidence Intervals (CIs) for the model coefficients.
- Calculate the p-values for hypothesis testing.
- Exponentiate the coefficients to obtain the Relative Risk Ratios (RRRs).
3.5.1 Calculate 95% Confidence Intervals and Combine with Coefficients
# Extract coefficients from the model
b_fitmlog2 <- coef(fitmlog2)
# Calculate confidence intervals
ci_fitmlog2 <- confint(fitmlog2)
# Combine coefficients and confidence intervals
b_ci_fitmlog2 <- cbind(b_fitmlog2, ci_fitmlog2)
b_ci_fitmlog2 b_fitmlog2 2.5 % 97.5 %
(Intercept):1 -1.1786550 -2.29930795 -0.05800204
(Intercept):2 -2.5649494 -4.59888160 -0.53101711
detc22:1 -0.3925617 -1.63588117 0.85075776
detc22:2 0.7060506 -1.41689336 2.82899454
detc23:1 0.1978257 -0.96565093 1.36130242
detc23:2 2.1059956 0.05519791 4.15679321-
coef(fitmlog2): Extracts the estimated coefficients from the modelfitmlog2. -
confint(fitmlog2): Calculates the 95% confidence intervals for the model coefficients. -
cbind(b_fitmlog2, ci_fitmlog2): Combines the coefficients and their confidence intervals into a single table. -
b_fitmlog2: The estimated coefficients. -
2.5 %and97.5 %: The lower and upper bounds of the 95% confidence intervals for each coefficient.
3.5.2 Exponentiate Coefficients to Obtain Relative Risk Ratios (RRRs)
# Exponentiate the combined coefficients and confidence intervals
rrr_fitmlog2 <- exp(b_ci_fitmlog2)
# Combine the RRRs with the original table
tab_fitmlog2 <- cbind(b_ci_fitmlog2, rrr_fitmlog2)
# Name the columns of the combined table
colnames(tab_fitmlog2) <- c("b", "lower b", "upper b", "rrr", "lower rrr", "upper rrr")
# Display the combined table
tab_fitmlog2 b lower b upper b rrr lower rrr
(Intercept):1 -1.1786550 -2.29930795 -0.05800204 0.30769231 0.10032825
(Intercept):2 -2.5649494 -4.59888160 -0.53101711 0.07692308 0.01006308
detc22:1 -0.3925617 -1.63588117 0.85075776 0.67532468 0.19478066
detc22:2 0.7060506 -1.41689336 2.82899454 2.02597403 0.24246610
detc23:1 0.1978257 -0.96565093 1.36130242 1.21875000 0.38073529
detc23:2 2.1059956 0.05519791 4.15679321 8.21527778 1.05674974
upper rrr
(Intercept):1 0.9436480
(Intercept):2 0.5880066
detc22:1 2.3414204
detc22:2 16.9284313
detc23:1 3.9012711
detc23:2 63.8663880-
exp(b_ci_fitmlog2): Exponentiates the coefficients and their confidence intervals to obtain the RRRs. -
cbind(b_ci_fitmlog2, rrr_fitmlog2): Combines the original table with the RRRs. -
colnames(tab_fitmlog2) <- c('b', 'lower b', 'upper b', 'rrr', 'lower rrr', 'upper rrr'): Renames the columns of the combined table for clarity. -
b: The estimated coefficients. -
lower bandupper b: The lower and upper bounds of the 95% confidence intervals for the coefficients. -
rrr: The Relative Risk Ratios (exponentiated coefficients). -
lower rrrandupper rrr: The lower and upper bounds of the 95% confidence intervals for the RRRs.
Interpretation
-
(Intercept):1
- The log-odds of being in the “over a year ago” category relative to “never” are -1.1786550 when all predictors are at their reference levels.
- The confidence interval for this log-odds is [-2.29930795, -0.0580204], indicating it is statistically significant since it does not include zero.
- The relative risk ratio (RRR) is 0.3076921, meaning that the likelihood of being in the “over a year ago” category is about 0.31 times that of the “never” category.
- The 95% confidence interval for the RRR is [0.10032825, 0.9434680], indicating the precision of the estimate.
-
detc22:1
- The log-odds of being in the “over a year ago” category relative to “never” decrease by 0.3925617 for each one-unit increase in
detc22. - The confidence interval for this coefficient is [-1.63588117, 0.85075776], which includes zero, indicating it is not statistically significant.
- The RRR is 0.67532468, meaning the likelihood decreases to about 0.68 times for each one-unit increase in
detc22. - The 95% confidence interval for the RRR is [0.19478066, 2.3414204].
- The log-odds of being in the “over a year ago” category relative to “never” decrease by 0.3925617 for each one-unit increase in
-
detc23:1
- The log-odds of being in the “over a year ago” category relative to “never” increase by 0.1978257 for each one-unit increase in
detc23. - The confidence interval for this coefficient is [-0.96556939, 1.36130242], which includes zero, indicating it is not statistically significant.
- The RRR is 1.21875000, meaning the likelihood increases to about 1.22 times for each one-unit increase in
detc23. - The 95% confidence interval for the RRR is [0.38073529, 3.9012711].
- The log-odds of being in the “over a year ago” category relative to “never” increase by 0.1978257 for each one-unit increase in
-
(Intercept):2
- The log-odds of being in the “within a year” category relative to “never” are -2.5649494 when all predictors are at their reference levels.
- The confidence interval for this log-odds is [-4.59888160, -0.5310171], indicating it is statistically significant.
- The RRR is 0.07692308, meaning that the likelihood of being in the “within a year” category is about 0.08 times that of the “never” category.
- The 95% confidence interval for the RRR is [0.01006308, 0.5888066].
-
detc22:2
- The log-odds of being in the “within a year” category relative to “never” increase by 0.7060506 for each one-unit increase in
detc22. - The confidence interval for this coefficient is [-1.41689336, 2.82899454], which includes zero, indicating it is not statistically significant.
- The RRR is 2.02597403, meaning the likelihood increases to about 2.03 times for each one-unit increase in
detc22. - The 95% confidence interval for the RRR is [0.24246610, 16.9284133].
- The log-odds of being in the “within a year” category relative to “never” increase by 0.7060506 for each one-unit increase in
-
detc23:2
- The log-odds of being in the “within a year” category relative to “never” increase by 2.1059956 for each one-unit increase in
detc23. - The confidence interval for this coefficient is [0.05519791, 4.15679321], indicating it is statistically significant since it does not include zero.
- The RRR is 8.21527778, meaning the likelihood increases to about 8.22 times for each one-unit increase in
detc23. - The 95% confidence interval for the RRR is [1.05674974, 63.8663880], indicating a high degree of variability in the estimate.
- The log-odds of being in the “within a year” category relative to “never” increase by 2.1059956 for each one-unit increase in
4 Practical 2: Additional Covariates
In this exercise, the objective is to extend the initial multinomial logistic regression model by adding more covariates to better understand the relationships between the predictors and the outcome variable. This helps in refining the model and potentially improving its explanatory power.
4.1 Estimation
Define the additional covariates to be included in the model. These include variables such as:
-
sympt(perception of symptoms) -
pb(perceived benefit) -
hist(family history) -
bse(breast self-exam training) -
detc(detection likelihood)
4.1.1 Defining the Model
The model is specified as a multinomial logistic regression, with the dataset dat.m3. The additional covariates are included in the model.
# Fitting the model with additional covariates
fitmlog3 <- vglm(me3 ~ factor(sympt) + pb + hist2 + bse + detc2,
family = multinomial(),
data = dat.m3
)
summary(fitmlog3)
Call:
vglm(formula = me3 ~ factor(sympt) + pb + hist2 + bse + detc2,
family = multinomial(), data = dat.m3)
Coefficients:
Estimate Std. Error z value Pr(>|z|)
(Intercept):1 -0.98609 1.11184 -0.887 0.37513
(Intercept):2 -2.99875 1.53905 NA NA
factor(sympt)2:1 -0.29008 0.64406 -0.450 0.65243
factor(sympt)2:2 0.11004 0.92273 0.119 0.90508
factor(sympt)3:1 0.81731 0.53979 1.514 0.12999
factor(sympt)3:2 1.92471 0.77757 2.475 0.01331 *
factor(sympt)4:1 1.13224 0.54767 2.067 0.03870 *
factor(sympt)4:2 2.45699 0.77530 3.169 0.00153 **
pb:1 -0.14821 0.07637 -1.941 0.05230 .
pb:2 -0.21944 0.07551 -2.906 0.00366 **
hist2yes:1 1.06544 0.45940 2.319 0.02038 *
hist2yes:2 1.36624 0.43752 3.123 0.00179 **
bse:1 1.05214 0.51499 2.043 0.04105 *
bse:2 1.29167 0.52988 2.438 0.01478 *
detc22:1 -0.92439 0.71375 -1.295 0.19528
detc22:2 0.01702 1.16169 0.015 0.98831
detc23:1 -0.69053 0.68712 -1.005 0.31491
detc23:2 0.90414 1.12661 0.803 0.42225
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Names of linear predictors: log(mu[,1]/mu[,3]), log(mu[,2]/mu[,3])
Residual deviance: 693.9019 on 806 degrees of freedom
Log-likelihood: -346.951 on 806 degrees of freedom
Number of Fisher scoring iterations: 5
Warning: Hauck-Donner effect detected in the following estimate(s):
'(Intercept):2'
Reference group is level 3 of the response-
vglm(): Fits a multinomial logistic regression model withme3as the dependent variable and the specified covariates as predictors. -
factor(sympt)convertssymptto a factor variable.
4.1.2 Checking and Recoding sympt Variable
The variable sympt has 4 categories, hence it is needed to be recoded into a binary factor variable symptd for simplicity.
# A tibble: 4 × 2
sympt n
<dbl> <int>
1 1 40
2 2 73
3 3 160
4 4 139# recode `sympt` into a binary variable `symptd`
dat.m3 <- dat.m3 %>%
mutate(symptd = ifelse(sympt < 3, "<3", "3 or more"))
dat.m3 %>% count(sympt, symptd)# A tibble: 4 × 3
sympt symptd n
<dbl> <chr> <int>
1 1 <3 40
2 2 <3 73
3 3 3 or more 160
4 4 3 or more 139-
mutate(symptd = ifelse(sympt < 3, '<3', '3 or more')): Recodesymptinto a binary variablesymptd, where categories less than 3 are labeled<3and categories 3 or more are labeled3 or more.
4.1.3 Refit the Model with New symptd Variable
# Refit the model using the new binary variable `symptd`
fitmlog4 <- vglm(me3 ~ symptd + pb + hist2 + bse + detc2,
family = multinomial(),
data = dat.m3
)
summary(fitmlog4)
Call:
vglm(formula = me3 ~ symptd + pb + hist2 + bse + detc2, family = multinomial(),
data = dat.m3)
Coefficients:
Estimate Std. Error z value Pr(>|z|)
(Intercept):1 -0.99877 1.07197 -0.932 0.351485
(Intercept):2 -2.70375 1.43422 -1.885 0.059406 .
symptd3 or more:1 1.12136 0.35720 3.139 0.001693 **
symptd3 or more:2 2.09534 0.45739 4.581 4.62e-06 ***
pb:1 -0.16811 0.07417 -2.266 0.023426 *
pb:2 -0.25101 0.07293 -3.442 0.000578 ***
hist2yes:1 1.01406 0.45381 2.235 0.025446 *
hist2yes:2 1.29328 0.43353 2.983 0.002853 **
bse:1 1.02859 0.51398 2.001 0.045366 *
bse:2 1.24397 0.52630 2.364 0.018097 *
detc22:1 -0.90213 0.71463 -1.262 0.206813
detc22:2 0.09028 1.16079 0.078 0.938011
detc23:1 -0.66982 0.68759 -0.974 0.329979
detc23:2 0.97281 1.12604 0.864 0.387627
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Names of linear predictors: log(mu[,1]/mu[,3]), log(mu[,2]/mu[,3])
Residual deviance: 697.4959 on 810 degrees of freedom
Log-likelihood: -348.748 on 810 degrees of freedom
Number of Fisher scoring iterations: 5
No Hauck-Donner effect found in any of the estimates
Reference group is level 3 of the responseSummary of Output
-
Intercepts
- (Intercept):1: The log-odds of being in the “over a year ago” category relative to “never” are -0.99877 when all predictors are at their reference levels. This is not statistically significant (p = 0.351485).
- (Intercept):2: The log-odds of being in the “within a year” category relative to “never” are -2.70375 when all predictors are at their reference levels. This is marginally significant (p = 0.059406).
-
Predictors
symptd3 or more:1: The log-odds of being in the “over a year ago” category relative to “never” increase by 1.12136 for those with symptoms 3 or more compared to less than 3. This is statistically significant (p = 0.001693).
symptd3 or more:2: The log-odds of being in the “within a year” category relative to “never” increase by 2.09534 for those with symptoms 3 or more compared to less than 3. This is highly significant (p < 0.001).
pb:1: The log-odds of being in the “over a year ago” category relative to “never” decrease by 0.16811 for each unit increase in perceived benefit. This is statistically significant (p = 0.023426).
pb:2: The log-odds of being in the “within a year” category relative to “never” decrease by 0.25101 for each unit increase in perceived benefit. This is highly significant (p < 0.001).
hist2yes:1: The log-odds of being in the “over a year ago” category relative to “never” increase by 1.01406 for those with a family history of breast cancer compared to those without. This is statistically significant (p = 0.025446).
hist2yes:2: The log-odds of being in the “within a year” category relative to “never” increase by 1.29383 for those with a family history of breast cancer compared to those without. This is highly significant (p = 0.002835).
bse:1: The log-odds of being in the “over a year ago” category relative to “never” increase by 1.28509 for those with breast self-exam training compared to those without. This is statistically significant (p = 0.012473).
bse:2: The log-odds of being in the “within a year” category relative to “never” increase by 1.24397 for those with breast self-exam training compared to those without. This is statistically significant (p = 0.018097).
detc22:1: The log-odds of being in the “over a year ago” category relative to “never” decrease by 0.90213 for those in the second level of detection likelihood compared to the reference level. This is not statistically significant (p = 0.206813).
detc22:2: The log-odds of being in the “within a year” category relative to “never” increase by 0.09028 for those in the second level of detection likelihood compared to the reference level. This is not statistically significant (p = 0.938011).
detc23:1: The log-odds of being in the “over a year ago” category relative to “never” decrease by 0.66982 for those in the third level of detection likelihood compared to the reference level. This is not statistically significant (p = 0.329979).
detc23:2: The log-odds of being in the “within a year” category relative to “never” increase by 0.97281 for those in the third level of detection likelihood compared to the reference level. This is not statistically significant (p = 0.387627).
-
Interpretation
- The model includes additional covariates such as symptoms, perceived benefit, family history, breast self-exam training, and detection likelihood.
- Significant predictors include
symptd,pb,hist2, andbse, indicating their impact on the likelihood of having a mammogram within a year or over a year ago relative to never. - The model fit statistics suggest that the model adequately fits the data.
- Confidence intervals and RRRs provide further insights into the effects of these predictors.
4.1.4 Convert detc to a Binary Variable
The variable detc initially has 3 categories. The aim is to recode it into a binary variable detcd.
-
ifelse(detc < 3, '<3', '3 or more'): This function checks the value ofdetcfor each observation. - If
detcis less than 3, it assigns the value ‘<3’ todetcd. - If
detcis 3 or more, it assigns the value ‘3 or more’ todetcd. - This effectively recodes
detcfrom 3 categories into a binary variable with 2 categories: ‘<3’ and ‘3 or more’.
4.1.5 Refit the Model with the New detcd Variable
fitmlog5 <- vglm(me3 ~ symptd + pb + hist2 + bse + detcd,
family = multinomial(),
data = dat.m3
)
summary(fitmlog5)
Call:
vglm(formula = me3 ~ symptd + pb + hist2 + bse + detcd, family = multinomial(),
data = dat.m3)
Coefficients:
Estimate Std. Error z value Pr(>|z|)
(Intercept):1 -1.82388 0.85509 -2.133 0.032928 *
(Intercept):2 -2.62376 0.92639 -2.832 0.004622 **
symptd3 or more:1 1.12742 0.35636 3.164 0.001558 **
symptd3 or more:2 2.09475 0.45742 4.579 4.66e-06 ***
pb:1 -0.15432 0.07262 -2.125 0.033587 *
pb:2 -0.24947 0.07258 -3.437 0.000588 ***
hist2yes:1 1.06318 0.45284 2.348 0.018885 *
hist2yes:2 1.30986 0.43360 3.021 0.002520 **
bse:1 0.95601 0.50734 1.884 0.059515 .
bse:2 1.23701 0.52542 2.354 0.018557 *
detcd3 or more:1 0.11416 0.31821 0.359 0.719786
detcd3 or more:2 0.88518 0.35624 2.485 0.012962 *
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Names of linear predictors: log(mu[,1]/mu[,3]), log(mu[,2]/mu[,3])
Residual deviance: 699.1326 on 812 degrees of freedom
Log-likelihood: -349.5663 on 812 degrees of freedom
Number of Fisher scoring iterations: 5
No Hauck-Donner effect found in any of the estimates
Reference group is level 3 of the responseSummary of Output
-
Intercepts
- (Intercept):1: The log-odds of being in the “over a year ago” category relative to “never” are -1.82388 when all predictors are at their reference levels. This is statistically significant (p = 0.032928).
- (Intercept):2: The log-odds of being in the “within a year” category relative to “never” are -2.62376 when all predictors are at their reference levels. This is highly significant (p = 0.004622).
-
Predictors
symptd3 or more:1: The log-odds of being in the “over a year ago” category relative to “never” increase by 1.12742 for those with symptoms 3 or more compared to less than 3. This is statistically significant (p = 0.001558).
symptd3 or more:2: The log-odds of being in the “within a year” category relative to “never” increase by 2.09475 for those with symptoms 3 or more compared to less than 3. This is highly significant (p < 0.001).
pb:1: The log-odds of being in the “over a year ago” category relative to “never” decrease by 0.15432 for each unit increase in perceived benefit. This is statistically significant (p = 0.033587).
pb:2: The log-odds of being in the “within a year” category relative to “never” decrease by 0.24947 for each unit increase in perceived benefit. This is highly significant (p < 0.001).
hist2yes:1: The log-odds of being in the “over a year ago” category relative to “never” increase by 1.06318 for those with a family history of breast cancer compared to those without. This is statistically significant (p = 0.018885).
hist2yes:2: The log-odds of being in the “within a year” category relative to “never” increase by 1.38968 for those with a family history of breast cancer compared to those without. This is highly significant (p = 0.001520).
bse:1: The log-odds of being in the “over a year ago” category relative to “never” increase by 0.95601 for those with breast self-exam training compared to those without. This is marginally significant (p = 0.059515).
bse:2: The log-odds of being in the “within a year” category relative to “never” increase by 1.20137 for those with breast self-exam training compared to those without. This is statistically significant (p = 0.018557).
detcd3 or more:1: The log-odds of being in the “over a year ago” category relative to “never” increase by 0.11416 for those with detection likelihood 3 or more compared to less than 3. This is not statistically significant (p = 0.719786).
detcd3 or more:2: The log-odds of being in the “within a year” category relative to “never” increase by 0.88518 for those with detection likelihood 3 or more compared to less than 3. This is statistically significant (p = 0.012962).
-
Interpretation
- The model includes covariates such as symptoms, perceived benefit, family history, breast self-exam training, and detection likelihood.
- Significant predictors include
symptd,pb,hist2, andbse, indicating their impact on the likelihood of having a mammogram within a year or over a year ago relative to never. - The model fit statistics suggest that the model adequately fits the data.
- Confidence intervals and RRRs provide further insights into the effects of these predictors.
4.2 Inferences
- Retrieve the log-odds and RRRs (Relative Risk Ratios).
- Calculate the 95% Confidence Intervals (CIs) for both the log-odds and the RRRs.
- Combine these results into a table.
- Present the results in a table
4.2.1 Retrieve Log-Odds and RRRs
The log-odds are the coefficients from the multinomial logistic regression model, and the RRRs are the exponentiated coefficients.
-
coef(fitmlog5): Extracts the coefficients (log-odds) from the modelfitmlog5. -
exp(coef(fitmlog5)): Exponentiates the coefficients to obtain the RRRs. -
cbind(coef(fitmlog5), exp(coef(fitmlog5))): Combines the log-odds and RRRs into one table.
4.2.2 Calculate 95% Confidence Intervals
-
confint(fitmlog5): Computes the 95% CIs for the log-odds. -
exp(confint(fitmlog5)): Computes the 95% CIs for the RRRs by exponentiating the CIs for the log-odds. -
cbind(confint(fitmlog5), exp(confint(fitmlog5))): Combines the CIs for log-odds and RRRs into one table.
4.2.3 Combine the Results into a Table
-
cbind(b_rrr_fitmlog5, ci_b_rrr_fitmlog5): Combines the log-odds, RRRs, and their CIs. -
colnames(res_fitmlog5) <- c('b', 'rrr', 'lower 95% b', 'upper 95% b', 'lower 95% rrr', 'upper 95% rrr'): Renames the columns of the combined table for clarity.
4.2.4 Display the Results
kable(res_fitmlog5, digits = 3) %>%
kable_styling(bootstrap_options = c("striped", "hover", "condensed", "responsive"))| b | rrr | lower 95% b | upper 95% b | lower 95% rrr | upper 95% rrr | |
|---|---|---|---|---|---|---|
| (Intercept):1 | -1.824 | 0.161 | -3.500 | -0.148 | 0.030 | 0.862 |
| (Intercept):2 | -2.624 | 0.073 | -4.439 | -0.808 | 0.012 | 0.446 |
| symptd3 or more:1 | 1.127 | 3.088 | 0.429 | 1.826 | 1.536 | 6.208 |
| symptd3 or more:2 | 2.095 | 8.123 | 1.198 | 2.991 | 3.314 | 19.911 |
| pb:1 | -0.154 | 0.857 | -0.297 | -0.012 | 0.743 | 0.988 |
| pb:2 | -0.249 | 0.779 | -0.392 | -0.107 | 0.676 | 0.898 |
| hist2yes:1 | 1.063 | 2.896 | 0.176 | 1.951 | 1.192 | 7.034 |
| hist2yes:2 | 1.310 | 3.706 | 0.460 | 2.160 | 1.584 | 8.669 |
| bse:1 | 0.956 | 2.601 | -0.038 | 1.950 | 0.962 | 7.031 |
| bse:2 | 1.237 | 3.445 | 0.207 | 2.267 | 1.230 | 9.649 |
| detcd3 or more:1 | 0.114 | 1.121 | -0.510 | 0.738 | 0.601 | 2.091 |
| detcd3 or more:2 | 0.885 | 2.423 | 0.187 | 1.583 | 1.206 | 4.871 |
-
library(knitr)andlibrary(kableExtra): Load the necessary libraries for table formatting. -
kable(res_fitmlog5, digits = 3): Creates a table withres_fitmlog5data, formatting the values to three decimal places. -
kable_styling(bootstrap_options = c("striped", "hover", "condensed", "responsive")): Applies styling options to the table for better presentation.
Interpretation
-
(Intercept):1
- The log-odds of being in the “over a year ago” category relative to “never” are -1.824 when all predictors are at their reference levels.
- The RRR indicates that the odds of being in the “over a year ago” category are 0.161 times those of “never”.
- The CIs suggest this estimate is statistically significant.
-
symptd3 or more:1
- The log-odds of being in the “over a year ago” category relative to “never” increase by 1.127 for those with symptoms 3 or more compared to less than 3.
- The RRR indicates that the odds are 3.088 times higher.
- The CIs suggest this estimate is statistically significant.
-
pb:1
- The log-odds of being in the “over a year ago” category relative to “never” decrease by 0.154 for each unit increase in perceived benefit.
- The RRR indicates that the odds are 0.857 times lower.
- The CIs suggest this estimate is statistically significant.
-
hist2yes:1
- The log-odds of being in the “over a year ago” category relative to “never” increase by 1.063 for those with a family history of breast cancer compared to those without.
- The RRR indicates that the odds are 2.896 times higher.
- The CIs suggest this estimate is statistically significant.
-
bse:1
- The log-odds of being in the “over a year ago” category relative to “never” increase by 0.956 for those with breast self-exam training compared to those without.
- The RRR indicates that the odds are 2.601 times higher.
- The CIs suggest this estimate is marginally significant.
-
detcd3 or more:1
- The log-odds of being in the “over a year ago” category relative to “never” increase by 0.114 for those with detection likelihood 3 or more compared to less than 3.
- The RRR indicates that the odds are 1.121 times higher.
- The CIs suggest this estimate is not statistically significant.
-
(Intercept):2
- The log-odds of being in the “within a year” category relative to “never” are -2.624 when all predictors are at their reference levels.
- The RRR indicates that the odds of being in the “within a year” category are 0.073 times those of “never”.
- The CIs suggest this estimate is statistically significant.
-
symptd3 or more:2
- The log-odds of being in the “within a year” category relative to “never” increase by 2.095 for those with symptoms 3 or more compared to less than 3.
- The RRR indicates that the odds are 8.123 times higher.
- The CIs suggest this estimate is highly significant.
-
pb:2
- The log-odds of being in the “within a year” category relative to “never” decrease by 0.249 for each unit increase in perceived benefit.
- The RRR indicates that the odds are 0.779 times lower.
- The CIs suggest this estimate is highly significant.
-
hist2yes:2
- The log-odds of being in the “within a year” category relative to “never” increase by 1.310 for those with a family history of breast cancer compared to those without.
- The RRR indicates that the odds are 3.706 times higher.
- The CIs suggest this estimate is statistically significant.
-
bse:2
- The log-odds of being in the “within a year” category relative to “never” increase by 1.237 for those with breast self-exam training compared to those without.
- The RRR indicates that the odds are 3.445 times higher.
- The CIs suggest this estimate is statistically significant.
-
detcd3 or more:2
- The log-odds of being in the “within a year” category relative to “never” increase by 0.885 for those with detection likelihood 3 or more compared to less than 3.
- The RRR indicates that the odds are 2.423 times higher.
- The CIs suggest this estimate is statistically significant.
4.3 Prediction
4.3.1 Predict the Log Odds
Predict the log odds for the given multinomial logistic regression model (fitmlog1).
# Predicts the log odds and display the first 6 observations
head(predict.vgam(fitmlog1, type = 'link')) log(mu[,1]/mu[,3]) log(mu[,2]/mu[,3])
1 -1.2504928 -0.9509763
2 -1.2504928 -0.9509763
3 -0.2411621 0.3053816
4 -1.2504928 -0.9509763
5 -1.2504928 -0.9509763
6 -1.2504928 -0.9509763-
predict.vgam(fitmlog1, type = 'link'): Predicts the log odds (logit values) for each observation using the modelfitmlog1. Thetype = 'link'option specifies that the predictions should be on the logit (log-odds) scale. -
head(): Displays the first 6 observations. - log(mu[,1]/mu[,3]): The predicted log odds of being in the “over a year ago” category vs. “never”.
- log(mu[,2]/mu[,3]): The predicted log odds of being in the “within a year” category vs. “never”.
4.3.1.1 Manual Calculation of Log Odds
Manually verify these predictions using the coefficients from the model and the values of the predictors for that observation.
- (Intercept):1: -1.2504928
- (Intercept):2: -0.9509763
- hist2yes:1: 1.0093308
- hist2yes:2: 1.2565379
head(dat.m1[1:3,])# A tibble: 3 × 7
obs me sympt pb hist bse detc
<dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
1 1 0 3 7 0 1 2
2 2 0 2 11 0 1 3
3 3 0 3 8 1 1 34.3.1.1.1 1st Observation (hist2 = no)
Logit for “over a year ago” ([1]) vs “never” ([3]):
\[ \text{logit}(P(\text{me3} = 1)) = \beta_0^{[1]} + \beta_{\text{hist2yes}}^{[1]} \cdot \text{hist2} \] \[ = -1.2504928 + 1.0093308 \cdot 0 \] \[ = -1.2504928 + 0 \] \[ = -1.2504928 \]
Logit for “within a year” ([2]) vs “never” ([3]):
\[ \text{logit}(P(\text{me3} = 2)) = \beta_0^{[2]} + \beta_{\text{hist2yes}}^{[2]} \cdot \text{hist2} \] \[ = -0.9509763 + 1.2565379 \cdot 0 \] \[ = -0.9509763 + 0 \] \[ = -0.9509763 \]
4.3.1.1.2 3rd Observation (hist2 = yes)
Logit for “over a year ago” ([1]) vs “never” ([3]):
\[ \text{logit}(P(\text{me3} = 1)) = \beta_0^{[1]} + \beta_{\text{hist2yes}}^{[1]} \cdot \text{hist2} \] \[ = -1.2504928 + 1.0093308 \cdot 1 \] \[ = -1.2504928 + 1.0093308 \] \[ = -0.241162 \]
Logit for “within a year” ([2]) vs “never” ([3]):
\[ \text{logit}(P(\text{me3} = 2)) = \beta_0^{[2]} + \beta_{\text{hist2yes}}^{[2]} \cdot \text{hist2} \] \[ = -0.9509763 + 1.2565379 \cdot 1 \] \[ = -0.9509763 + 1.2565379 \] \[ = 0.3053816 \]
4.3.2 Predict the Probability
Using fitmlog1 model, the log odds for the 1st observation:
- Log odds for “over a year ago”: -1.2504928
- Log odds for “within a year”: -0.9509763
# Predicts the probabilities for the first 6 observations
head(predict.vgam(fitmlog1, type = 'response')) over.a.year.ago within.a.year never
1 0.1711957 0.2309783 0.5978261
2 0.1711957 0.2309783 0.5978261
3 0.2500000 0.4318182 0.3181818
4 0.1711957 0.2309783 0.5978261
5 0.1711957 0.2309783 0.5978261
6 0.1711957 0.2309783 0.5978261-
predict.vgam(fitmlog1, type = 'response'): Predicts the probabilities for each outcome category for the observations using the modelfitmlog1. Thetype = 'response'option specifies that the predictions should be on the probability scale. -
head(): Displays the first 6 observations. - over.a.year.ago: Probability of being in the “over a year ago” category.
- within.a.year: Probability of being in the “within a year” category.
- never: Probability of being in the “never” category.
The calculations below demonstrate how to use the predicted log odds to derive the probabilities for each outcome category, ensuring that the probabilities sum to 1 for each observation.
4.3.2.1 Probability of Being in the Reference Group (“Never”)
\[ P(\text{me3} = \text{never}) = \frac{1}{1 + \exp(\text{logit}_{\text{over a year ago}}) + \exp(\text{logit}_{\text{within a year}})} \] \[ P(\text{me3} = \text{never}) = \frac{1}{1 + \exp(-1.2504928) + \exp(-0.9509763)} \] \[ = \frac{1}{1 + 0.2865693 + 0.3867897} \] \[ = \frac{1}{1.673359} \] \[ = 0.5978261 \]
4.3.2.2 Probability of Being in the “Over a Year Ago” Group
\[ P(\text{me3} = \text{over a year ago}) = \frac{\exp(\text{logit}_{\text{over a year ago}})}{1 + \exp(\text{logit}_{\text{over a year ago}}) + \exp(\text{logit}_{\text{within a year}})} \] \[ P(\text{me3} = \text{over a year ago}) = \frac{\exp(-1.2504928)}{1 + \exp(-1.2504928) + \exp(-0.9509763)} \] \[ = \frac{0.2865693}{1 + 0.2865693 + 0.3867897} \] \[ = \frac{0.2865693}{1.673359} \] \[ = 0.1711957 \]
4.3.2.3 Probability of Being in the “Within a Year” Group
\[ P(\text{me3} = \text{within a year}) = \frac{\exp(\text{logit}_{\text{within a year}})}{1 + \exp(\text{logit}_{\text{over a year ago}}) + \exp(\text{logit}_{\text{within a year}})} \] \[ P(\text{me3} = \text{within a year}) = \frac{\exp(-0.9509763)}{1 + \exp(-1.2504928) + \exp(-0.9509763)} \] \[ = \frac{0.3867897}{1 + 0.2865693 + 0.3867897} \] \[ = \frac{0.3867897}{1.673359} \] \[ = 0.2309783 \]
5 Practical 3: Alternative Approach
This exercise demonstrates an alternative approach to fitting a multinomial logistic regression model using the multinom function from the nnet package. It includes steps to fit the model, change the reference category, extract p-values, and compute confidence intervals.
In the VGAM::vglm function, the reference group is the category with the largest number of observations. In contrast, the nnet::multinom function uses the category with the smallest number of observations as the reference group.
5.1 Estimation
5.1.1 Fitting the Model
# Fit the multinomial logistic regression model
mlog_nnet1 <- multinom(me3 ~ hist2, data = dat.m3)# weights: 9 (4 variable)
initial value 452.628263
final value 396.169969
converged# Summary of the model
summary(mlog_nnet1)Call:
multinom(formula = me3 ~ hist2, data = dat.m3)
Coefficients:
(Intercept) hist2yes
within.a.year 0.2995173 0.2470274
never 1.2504941 -1.0093349
Std. Errors:
(Intercept) hist2yes
within.a.year 0.1662460 0.413737
never 0.1428933 0.427500
Residual Deviance: 792.3399
AIC: 800.3399 Note that the reference group here is “over.a.year.ago” instead of “never.”
5.1.2 Changing the Reference Category
# Create new dataset and change the reference category
dat.m4 <- dat.m3
dat.m4 <- dat.m4 %>% mutate(me4 = relevel(me3, ref = "never"))
# Recheck the order of categories
levels(dat.m4$me4)[1] "never" "over.a.year.ago" "within.a.year" # Fit the model with the new reference category
mlog_nnet2 <- multinom(me4~ hist2, data = dat.m4)# weights: 9 (4 variable)
initial value 452.628263
iter 10 value 396.169969
iter 10 value 396.169969
final value 396.169969
converged# Summary of the new model
summary(mlog_nnet2)Call:
multinom(formula = me4 ~ hist2, data = dat.m4)
Coefficients:
(Intercept) hist2yes
over.a.year.ago -1.2504803 1.009384
within.a.year -0.9509794 1.256531
Std. Errors:
(Intercept) hist2yes
over.a.year.ago 0.1428926 0.4275097
within.a.year 0.1277115 0.3746633
Residual Deviance: 792.3399
AIC: 800.3399 The reference category is changed to “never,” which affects the interpretation of the coefficients.
5.2 Inferences
5.2.1 Getting P-Values
# Calculate z-values
z.test <- summary(mlog_nnet2)$coefficients / summary(mlog_nnet2)$standard.errors
# Calculate p-values
p.val <- (1 - pnorm(abs(z.test), 0, 1)) * 2
colnames(p.val) <- c('p-val intercept', 'p-val hist')
p.val p-val intercept p-val hist
over.a.year.ago 0.000000e+00 0.0182218315
within.a.year 9.592327e-14 0.0007972127- z.test: z-values calculated by dividing coefficients by their standard errors.
- p.val: Two-tailed p-values calculated from the z-values.
5.2.2 Confidence Intervals
# Calculate confidence intervals
confint(mlog_nnet2, level = 0.95), , over.a.year.ago
2.5 % 97.5 %
(Intercept) -1.5305447 -0.9704159
hist2yes 0.1714807 1.8472880
, , within.a.year
2.5 % 97.5 %
(Intercept) -1.2012893 -0.7006696
hist2yes 0.5222045 1.9908576- Confidence intervals: Provide a range within which the true coefficient values are expected to lie with 95% confidence.
6 Acknowledgment
I would like to express my gratitude to Professor Dr. Kamarul Imran Musa, Medical Epidemiologist and Statistician, and Professor in Epidemiology and Biostatistics at the School of Medical Sciences, Universiti Sains Malaysia, for his teaching and guidance on the subject matter.
7 References
- Lecture on Multinomial Logistic Regression by Professor Dr. Kamarul Imran Musa.
- Applied Logistic Regression, Second Edition by Hosmer and Lemeshow, Chapter 8: Special Topics. UCLA Institute for Digital Research & Education.
- STAT 504: Analysis of Discrete Data, Lesson 11: Multinomial Logistic Regression. Penn State Eberly College of Science.
- R Data Analysis Examples: Multinomial Logistic Regression. UCLA Institute for Digital Research & Education.
- The mlogit Package: Multinomial Logit Models in R. Croissant, Y. (2021).
- Long, J. S., & Freese, J. “Models for Nominal Outcomes.” In Regression Models for Categorical Dependent Variables Using Stata (3rd ed.). Stata Press.
- Hosmer, D. W. Jr., Lemeshow, S., & Sturdivant, R. X. (2013). Special Topics. In Applied Logistic Regression (3rd ed.). Wiley.
8 R Codes
The R codes used in this analysis are available at the following GitHub repository: DrPH-Epidemiology-Revision. This repository includes all scripts and data necessary to replicate the analyses presented in this work.